% read files and create ETData Objects
if ~exist('e13', 'var')
    % change workfolder
    cd('D:\\Workspace\\Compensation\\Data\\CK50S');

    % create object
    e13 = ETData('ST20181113D');
    e15 = ETData('ST20181115D');
    e16 = ETData('ST20181116D');
    e19 = ETData('ST20181119D');
    e21 = ETData('ST20181121DC');

    % add compensation
    e21.comp([1 2 3]);

    f1 = e13.plot; f2 = e15.plot; f3 = e16.plot; f4 = e19.plot; f5 = e21.plot;
end

% model
if ~exist('B', 'var') || isempty(B)
    % sx = [1 3 4 5 11];
    % sx = [1 4 5 8 11];
    sx = [1 4 8 9 11];
    sy = 1;

    % [X1, Y1, U1, V1] = model(e13, [1, 14]);
    % [X2, Y2, U2, V2] = model(e15, [1, 18]);
    % [X3, Y3, U3, V3] = model(e16, [1, 17]);
    % [X4, Y4, U4, V4] = model(e19, [1, 17]);
    [X1, Y1, U1, V1] = model(e16, [1, 17]);
    [X2, Y2, U2, V2] = model(e19, [1, 17]);
    [X3, Y3, U3, V3] = model(e21, [1, 17]);
    % X1 = e13.Temp(:, sx); Y1 = e13.Eddy(:, 1);
    % X2 = e15.Temp(:, sx); Y2 = e15.Eddy(:, 1);
    % X3 = e16.Temp(:, sx); Y3 = e16.Eddy(:, 1);
    % X4 = e19.Temp(:, sx); Y4 = e19.Eddy(:, 1);

    [X, Y   ] = modelxy(X1, Y1, X2, Y2);
    [X, Y, B] = modelxy(X, Y, X3, Y3);
    % [X, Y, B] = modelxy(X, Y, X4, Y4);

    % [X, Y    ] = modelxy(X1(1:13, :), Y1(1:13, :), X2(1:17, :), Y2(1:17, :));
    % [X, Y,   ] = modelxy(X, Y, X3(1:16, :), Y3(1:16, :));
    % [X, Y, B1] = modelxy(X, Y, X4(1:16, :), Y4(1:16, :));

    % [X, Y    ] = modelxy(X1(14:end, :), Y1(14:end, :), X2(18:end, :), Y2(18:end, :));
    % [X, Y,   ] = modelxy(X, Y, X3(17:end, :), Y3(17:end, :));
    % [X, Y, B2] = modelxy(X, Y, X4(17:end, :), Y4(17:end, :));

    U = [U1; U2; U3];
    V = [V1; V2; V3];
    D = regress(V, U);

    B(1) = [];
end

% plot
if ~exist('l1', 'var')
    figure(f1); subplot(2, 1, 1);
    % x  = e13.Temp(:, sx);
    % c  = e13.Eddy(1, 1) - x(1, :) * B;
    % y  = x * B + c;
    % c  = e13.Eddy(14, 1) - x(14, :) * B2;
    % y2 = x(14:end, :) * B2 + c;
    % dt = (x(14, :) - x(13, :)) * (D - B);
    % y(14:end) = y(14:end) + dt;
    % disp(max(abs(e13.Eddy(:, 1) - y)));
    % l1 = plot(e13.Time, y, '--');

    figure(f2); subplot(2, 1, 1);
    % x  = e15.Temp(:, sx);
    % c  = e15.Eddy(1, 1) - x(1, :) * B;
    % y  = x * B + c;
    % c  = e15.Eddy(18, 1) - x(18, :) * B2;
    % y2 = x(18:end, :) * B2 + c;
    % dt = (x(18, :) - x(17, :)) * (D - B);
    % y(18:end) = y(18:end) + dt;
    % disp(max(abs(e15.Eddy(:, 1) - y)));
    % l2 = plot(e15.Time, y, '--');

    figure(f3); subplot(2, 1, 1);
    x  = e16.Temp(:, sx);
    c  = e16.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    % c  = e16.Eddy(17, 1) - x(17, :) * B2;
    % y2 = x(17:end, :) * B2 + c;
    dt = (x(17, :) - x(16, :)) * (D - B);
    y(17:end) = y(17:end) + dt;
    disp(max(abs(e16.Eddy(:, 1) - y)));
    l3 = plot(e16.Time, y, '--');

    figure(f4); subplot(2, 1, 1);
    x  = e19.Temp(:, sx);
    c  = e19.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    % c  = e19.Eddy(17, 1) - x(17, :) * B2;
    % y2 = x(17:end, :) * B2 + c;
    dt = (x(17, :) - x(16, :)) * (D - B);
    y(17:end) = y(17:end) + dt;
    disp(max(abs(e19.Eddy(:, 1) - y)));
    l4 = plot(e19.Time, y, '--');

    figure(f5); subplot(2, 1, 1);
    x  = e21.Temp(:, sx);
    c  = e21.Eddy(1, 1) - x(1, :) * B;
    y  = x * B + c;
    % c  = e21.Eddy(17, 1) - x(17, :) * B2;
    % y2 = x(17:end, :) * B2 + c;
    dt = (x(17, :) - x(16, :)) * (D - B);
    y(17:end) = y(17:end) + dt;
    disp(max(abs(e21.Eddy(:, 1) - y)));
    l5 = plot(e21.Time, y, '--');
end
